Figure 4 - IP-MS
Charlotte Soneson, Michael Stadler, Merle Skribbe
2024-06-21
Parameter values
params## $rds150
## [1] "data/ipms_150_sce.rds"
##
## $rds500
## [1] "data/ipms_500_sce.rds"
##
## $rdsdal2
## [1] "data/ipms_dal2_sce.rds"
##
## $complexes
## [1] "data/complexes.json"
##
## $baitclass
## [1] "data/ipms_bait_class.txt"
##
## $dnabinders
## [1] "data/Gene_list_DNAbinders_proteincoding.txt"
##
## $rnabinders
## [1] "data/Gene_list_RNAbinders_proteincoding.txt"
##
## $dbdtxt
## [1] "data/TF_DBD_annotation.txt"
##
## $idmap
## [1] "data/id_mapping_table.txt"Load required packages and helper functions
suppressPackageStartupMessages({
library(SingleCellExperiment)
library(ggplot2)
library(cowplot)
library(einprot)
library(dplyr)
library(forcats)
library(jsonlite)
library(ComplexHeatmap)
library(colorspace)
library(ggrepel)
library(DT)
library(stringr)
library(tidyr)
library(kableExtra)
library(patchwork)
})
## Source scripts with required helper functions and settings
source("params/plot_settings.R")
source("params/get_testres_function.R")
source("params/ipms_params.R")
source("params/mapping_functions.R")
## Heatmap sizes
w_ipmsHm <- 84 # heatmap body width (mm)
h_ipmsHm <- 178 # total heatmap height (mm)
fs <- 8 # small font size
fl <- 10 # large font size
ft <- 12 # title font sizeRead data
We load two SingleCellExperiment objects, containing data and results from the low- and high-salt experiments, respectively, as well as one object with results from a separate run of Dal2.
sce150 <- readRDS(params$rds150)
sce500 <- readRDS(params$rds500)
scedal2 <- readRDS(params$rdsdal2)We then extract the full list of baits, as well as a classification according to the experiment(s) where a certain bait was pulled down.
idmap <- read.delim(params$idmap)
baitclass <- read.delim(params$baitclass)
DT::datatable(baitclass, rownames = FALSE)In the following set of IPs, we identified the bait as not pulled down.
(nobait_150 <- baitclass$Gene_name[baitclass$class == "nobait"])## [1] "Cha4" "Cuf2" "Atf21" "Cbf12" "Sre2" "SPBC1348.12" "Untagged"We also read a list of known complexes, which will be used for annotation purposes.
complexes <- jsonlite::read_json(params$complexes, simplifyVector = TRUE)
complexes$`MBF transcription complex`$color <- main_colors[5]
complexes$`CCAAT-binding factor complex`$color <- main_colors[1]
complexes$`Atf1-Pcr1`$color <- main_colors[3]
names(complexes)## [1] "NuA4" "SAGA" "19S proteasome - base"
## [4] "Clr6 HDAC - complex I'" "SREBP-SCAP" "SPBC16G5.16-SPBC530.08"
## [7] "Atf1-Pcr1" "chaperonin-containing T-complex" "CCAAT-binding factor complex"
## [10] "Dal2" "Cyc8(Ssn6)-Tup1" "Rad24-Rad25"
## [13] "Ste11-Cdc2-Cig2" "Pap1-Prr1" "MBF transcription complex"Finally we load a table classifying each TF to its DNA-binding domain family.
dbd <- read.delim(params$dbdtxt) |>
dplyr::mutate(Gene_name = .capitalize(Gene_name))
DT::datatable(dbd, extensions = "FixedColumns", rownames = FALSE,
filter = list(position = "top", clear = FALSE),
options = list(scrollX = TRUE,
fixedColumns = list(leftColumns = 1),
search = list(regex = FALSE, caseInsensitive = TRUE),
pageLength = 10))Helper functions
We define a number of additional helper functions for data extraction and summarization, which will be used below.
filterForHeatmap <- function(tstatmat, adjpmat, logfcmat, adjpThreshold,
logfcThreshold, maxNbrTargets) {
## Filter - require significant in at least one comparison
tstatfilt <- tstatmat
tstatfilt[adjpmat >= adjpThreshold] <- NA
tstatfilt[logfcmat <= logfcThreshold] <- NA
# tstatfilt[tstatfilt <= 0] <- NA
## Some experiments pull down a large number of proteins
# print(table(colSums(!is.na(tstatfilt))))
## -> remove proteins that are only being pulled down in these experiments
keep <- which(colSums(!is.na(tstatfilt)) <= maxNbrTargets)
tstatfilt <- tstatfilt[rowSums(!is.na(tstatfilt[, keep])) > 0, ]
tstatfilt[is.na(tstatfilt)] <- 0
tstatfilt <- tstatfilt[rowSums(tstatfilt) > 0, ]
## Subset the t-statistic matrix to only the retained proteins
stopifnot(colnames(tstatfilt) == colnames(tstatmat))
tstatmat <- tstatmat[rownames(tstatfilt), colnames(tstatfilt)]
colnames(tstatmat) <- .getSimplifiedComparison(colnames(tstatmat))
list(tstats = tstatmat, tstatsfilt = tstatfilt)
}
makeComplexAnnotation <- function(tstatmat, complexes, idmap,
show_legend = TRUE) {
rowannot <- data.frame(Complex = rep(NA, nrow(tstatmat)),
row.names = rownames(tstatmat))
rowcols <- list(Complex = c())
for (cmplx in names(complexes)) {
rowannot$Complex[.getPomBaseIdFromProtein(rownames(rowannot),
idmap = idmap) %in%
complexes[[cmplx]]$pombase_ids] <- cmplx
rowcols$Complex[cmplx] <- complexes[[cmplx]]$color
}
rowAnnotation(df = rowannot,
col = rowcols,
annotation_name_gp = gpar(fontsize = fl),
simple_anno_size_adjust = TRUE,
na_col = bg_color,
show_annotation_name = FALSE,
show_legend = show_legend,
annotation_legend_param = list(title_gp = gpar(fontsize = fl)))
}
makeComplexAndDBDAnnotation <- function(tstatmat, complexes, dbd, idmap) {
rowannot <- data.frame(Complex = rep(NA, nrow(tstatmat)),
row.names = rownames(tstatmat))
rowcols <- list(Complex = c())
for (cmplx in names(complexes)) {
rowannot$Complex[.getPomBaseIdFromProtein(rownames(rowannot),
idmap = idmap) %in%
complexes[[cmplx]]$pombase_ids] <- cmplx
rowcols$Complex[cmplx] <- complexes[[cmplx]]$color
}
rowannot$`DBD fam.` <- dbd$DBD_class[match(
.getPomBaseIdFromProtein(rownames(tstatmat), idmap = idmap),
dbd$PomBaseID)]
setfam <- c("KilA-N", "Zn(II)2Cys6", "bZIP")
unsetfam <- setdiff(unique(dbd$DBD_class), setfam)
rowannot$`DBD fam.`[rowannot$`DBD fam.` %in% unsetfam] <- "Other (combined)"
rowcols <- c(rowcols, list(`DBD fam.` = c(`KilA-N` = main_colors[5],
`Zn(II)2Cys6` = main_colors[4],
bZIP = main_colors[3],
`Other (combined)` = bg_color)))
rowAnnotation(df = rowannot,
col = rowcols,
annotation_name_gp = gpar(fontsize = fl),
simple_anno_size_adjust = TRUE,
na_col = bg_color,
show_annotation_name = TRUE,
show_legend = FALSE,
annotation_legend_param = list(title_gp = gpar(fontsize = fl)))
}
makeColLabels <- function(tstatmat, colLabel, fontsize) {
columnAnnotation(
TFnames = anno_mark(which = "column", side = "bottom",
at = match(colLabel, colnames(tstatmat)),
labels = .getProteinNameFromSimplifiedComparison(colLabel),
labels_gp = gpar(fontsize = fontsize)))
}
makeRowLabels <- function(tstatmat, rowLabel, fontsize) {
rowAnnotation(
TFnames = anno_mark(which = "row", side = "left",
at = match(rowLabel, rownames(tstatmat)),
labels = rowLabel,
labels_gp = gpar(fontsize = fontsize)))
}
makeHeatmapCol <- function(tstatmat, stringency = "high") {
if (stringency == "high") {
circlize::colorRamp2(
breaks = seq(3.0, 13.0, length.out = 64),
colors = colorRampPalette(binary_heatmap_colors[c("FALSE", "TRUE")])(64))
} else {
circlize::colorRamp2(
breaks = seq(1.5, 13.0, length.out = 64),
colors = colorRampPalette(binary_heatmap_colors[c("FALSE", "TRUE")])(64))
}
}
reorderColsByRows <- function(tstat) {
## Row names should be protein names, column names simplified comparison names
ord <- match(rownames(tstat), .getProteinNameFromSimplifiedComparison(colnames(tstat)))
ord <- ord[!is.na(ord)]
ord <- c(ord, setdiff(seq_len(ncol(tstat)), ord))
tstat[, ord]
}
makeHeatmapData <- function(sce, adjpthr, log2fcthr, conc) {
tc <- .getTestCols(sce, adjp_cutoff = adjpthr, logfc_cutoff = log2fcthr)
tstats <- filterForHeatmap(tstatmat = tc$tstat, adjpmat = tc$adjp,
logfcmat = tc$logfc, logfcThreshold = log2fcthr,
adjpThreshold = adjpthr, maxNbrTargets = Inf)
tstat <- tstats$tstats
## tstats$tstatsfile: protein x experiment
## Number of TFs that are pulled down by each bait
all_tfs <- setdiff(.getProteinNameFromComparison(colnames(tstats$tstatsfilt)),
c("Untagged", "untagged", NA))
nbrInt <- colSums(tstats$tstatsfilt[rownames(tstats$tstatsfilt) %in% all_tfs, ] > 0,
na.rm = TRUE)
names(nbrInt) <- .getSimplifiedComparison(names(nbrInt))
## Split by family - return a matrix with several columns:
## #pulled down TFs (including itself), #pulled down TFs (excluding itself),
## #pulled down TFs from the same family (excluding itself),
## #pulled down TFs from another family
tmptf <- tstats$tstatsfilt[match(.getProteinNameFromComparison(colnames(tstats$tstatsfilt)),
rownames(tstats$tstatsfilt)), ]
rownames(tmptf) <- .getProteinNameFromComparison(colnames(tstats$tstatsfilt))
dbdtf <- data.frame(PomBaseID = .getPomBaseIdFromComparison(colnames(tstats$tstatsfilt),
idmap = idmap),
bait = .getProteinNameFromComparison(colnames(tstats$tstatsfilt)),
DBD_class = NA)
for (i in seq_len(nrow(dbdtf))) {
lookfor <- dbdtf$PomBaseID[i]
if (!(lookfor %in% c("untagged", "Untagged"))) {
dbdtf$DBD_class[i] <- dbd$DBD_class[match(lookfor, dbd$PomBaseID)]
} else {
dbdtf$DBD_class[i] <- "no_class"
}
}
stopifnot(dbdtf$bait == .getProteinNameFromComparison(colnames(tmptf)))
nbrIntMat <- do.call(dplyr::bind_rows, lapply(seq_along(colnames(tmptf)), function(i) {
expr <- colnames(tmptf)[i]
fam <- dbdtf$DBD_class[i]
data.frame(experiment = expr,
bait = dbdtf$bait[i],
dbdfam = fam,
TFs = paste(rownames(tmptf)[-i][which(tmptf[-i, expr] > 0)],
collapse = ","),
nbrTFs = sum(tmptf[, expr] > 0, na.rm = TRUE),
nbrTFs_wobait = sum(tmptf[-i, expr] > 0, na.rm = TRUE),
nbrTFs_samefam = sum(tmptf[-i, expr] > 0 &
dbdtf$DBD_class[-i] == fam,
na.rm = TRUE),
nbrTFs_difffam = sum(tmptf[-i, expr] > 0 &
dbdtf$DBD_class[-i] != fam,
na.rm = TRUE)
)
}))
## Split experiments into three groups
if (conc == 150) {
stopifnot(all(colnames(tstat) == .getSimplifiedComparison(colnames(tstats$tstatsfilt))))
tstat_no_bait <-
tstat[, colSums(tstats$tstatsfilt) == 0 |
(.getOrigBaitNameFromComparison(colnames(tstats$tstatsfilt), idmap = idmap) %in%
baitclass$Gene_name[baitclass$class == "nobait"])]
tstat_150only <-
tstat[, colSums(tstats$tstatsfilt) > 0 &
(.getOrigBaitNameFromComparison(colnames(tstats$tstatsfilt), idmap = idmap) %in%
baitclass$Gene_name[baitclass$class == "150only"])]
tstat_150and500 <-
tstat[, colSums(tstats$tstatsfilt) > 0 &
(.getOrigBaitNameFromComparison(colnames(tstats$tstatsfilt), idmap = idmap) %in%
baitclass$Gene_name[baitclass$class == "150and500"])]
proteins_150only <- .getProteinNameFromSimplifiedComparison(colnames(tstat_150only))
stopifnot(all(proteins_150only %in% rownames(tstat)))
tstat_150and500 <- tstat_150and500[!(rownames(tstat_150and500) %in% proteins_150only), ]
## Cluster proteins based on the correlation between their t-stat vectors
dst <- sqrt(2 * (1 - cor(t(tstat_150and500), use = "pairwise.complete")))
dst[is.na(dst)] <- sqrt(2 * (1 - (-1)))
clst <- hclust(as.dist(dst))
## Reorder proteins by the clustering
protorder <- c(rownames(tstat_150and500)[clst$order], proteins_150only)
tstat <- tstat[protorder, ]
## Reorder experiments to keep a "diagonal" in the heatmap
tstat <- reorderColsByRows(tstat)
tmpbaits <- .getOrigBaitNameFromSimplifiedComparison(colnames(tstat),
idmap = idmap)
colsplit <- ifelse(
tmpbaits %in% baitclass$Gene_name[baitclass$class == "nobait"],
"No\nbait",
ifelse(
tmpbaits %in% baitclass$Gene_name[baitclass$class == "150only"],
"150 mM NaCl",
ifelse(
tmpbaits %in% baitclass$Gene_name[baitclass$class == "150and500"],
"150 and\n500 mM NaCl", "NA"
)
)
)
} else if (conc == 500) {
tstat_onlybait <-
tstat[, colSums(tstats$tstatsfilt > 0) == 1 |
(.getOrigBaitNameFromComparison(colnames(tstats$tstatsfilt),
idmap = idmap) %in%
baitclass$Gene_name[baitclass$class500 == "Lost in 500 mM"])]
tstat_several <-
tstat[, colSums(tstats$tstatsfilt > 0) > 1 &
(.getOrigBaitNameFromComparison(colnames(tstats$tstatsfilt),
idmap = idmap) %in%
baitclass$Gene_name[baitclass$class500 != "Lost in 500 mM"])]
tstat_onlybait <- tstat_onlybait[, colnames(tstat_onlybait) != "Untagged_plate"]
tstat_several <- tstat_several[, colnames(tstat_several) != "Untagged_plate"]
proteins_onlybait <-
.getProteinNameFromSimplifiedComparison(colnames(tstat_onlybait))
stopifnot(all(proteins_onlybait %in% rownames(tstat)))
tstat_several <- tstat_several[!(rownames(tstat_several) %in% proteins_onlybait), ]
## Cluster proteins based on the correlation between their t-stat vectors
dst <- sqrt(2 * (1 - cor(t(tstat_several), use = "pairwise.complete")))
dst[is.na(dst)] <- sqrt(2 * (1 - (-1)))
clst <- hclust(as.dist(dst))
## Reorder proteins by the clustering
protorder <- c(rownames(tstat_several)[clst$order], proteins_onlybait)
tstat <- tstat[protorder, ]
## Reorder conditions to keep a "diagonal" in the heatmap
tstat <- reorderColsByRows(tstat)
colsplit <- factor(
structure(baitclass$class500[
match(.getOrigBaitNameFromSimplifiedComparison(colnames(tstat),
idmap = idmap),
baitclass$Gene_name)],
names = colnames(tstat)),
levels = c("Retained in 500 mM", "Lost in 500 mM", ""))
levels(colsplit) <- c("Interactions\nin 500 mM NaCl",
"No interactions\nin 500 mM NaCl", "")
} else {
stop("Unknown concentration")
}
return(list(tstat = tstat, colsplit = colsplit, nbrIntMat = nbrIntMat))
}
makeBaitHeatmapData <- function(sce, adjpthr, log2fcthr, nbrIntMat) {
tc <- .getTestCols(sce, adjp_cutoff = adjpthr, logfc_cutoff = log2fcthr)$tstat
colnames(tc) <- .getSimplifiedComparison(colnames(tc))
shared_names <- intersect(rownames(tc),
.getProteinNameFromSimplifiedComparison(colnames(tc)))
shared_names <- setdiff(shared_names, .getProteinFromOrigBait(nobait_150, idmap = idmap))
stopifnot(all(shared_names %in% nbrIntMat$bait))
nbrIntMat$complex <- NA
for (cplx in c("MBF transcription complex", "CCAAT-binding factor complex", "Atf1-Pcr1")) {
nbrIntMat$complex[.getPomBaseIdFromProtein(nbrIntMat$bait,
idmap = idmap) %in%
complexes[[cplx]]$pombase_ids] <- cplx
}
shared_names <- nbrIntMat[match(shared_names, nbrIntMat$bait), ] |>
arrange(dbdfam, complex) |>
pull(bait)
mat <- tc[shared_names, match(shared_names, .getProteinNameFromSimplifiedComparison(colnames(tc)))]
rownames(nbrIntMat) <- .getSimplifiedComparison(nbrIntMat$experiment)
list(mat = mat, nbrIntMat = nbrIntMat)
}Extract data for heatmaps
hmdata_150 <- makeHeatmapData(sce = sce150, adjpthr = adjpThreshold,
log2fcthr = log2fcThreshold, conc = 150)
hmdata_150_loose <- makeHeatmapData(sce = sce150, adjpthr = adjpThresholdLoose,
log2fcthr = log2fcThreshold, conc = 150)
hmdata_500 <- makeHeatmapData(sce = sce500, adjpthr = adjpThreshold,
log2fcthr = log2fcThreshold, conc = 500)
hmdata_500_loose <- makeHeatmapData(sce = sce500, adjpthr = adjpThresholdLoose,
log2fcthr = log2fcThreshold, conc = 500)
baitdata_150 <- makeBaitHeatmapData(sce = sce150, adjpthr = adjpThreshold,
log2fcthr = log2fcThreshold,
nbrIntMat = hmdata_150$nbrIntMat)
baitdata_150_loose <- makeBaitHeatmapData(sce = sce150, adjpthr = adjpThresholdLoose,
log2fcthr = log2fcThreshold,
nbrIntMat = hmdata_150_loose$nbrIntMat)
baitdata_500 <- makeBaitHeatmapData(sce = sce500, adjpthr = adjpThreshold,
log2fcthr = log2fcThreshold,
nbrIntMat = hmdata_500$nbrIntMat)
baitdata_500_loose <- makeBaitHeatmapData(sce = sce500, adjpthr = adjpThresholdLoose,
log2fcthr = log2fcThreshold,
nbrIntMat = hmdata_500_loose$nbrIntMat)Export data for TFexplorer
Here, we export data for use with TFexplorer, for both the low- and high-salt conditions.
## Heatmap data
.prepDataForTFexplorer <- function(sce, outfileHeatmap, outfileVolcano, idmap) {
## Heatmap
## -------------------------------------------------------------------------
testres <- .getTestCols(sce = sce)
testres$tstat <- testres$tstat[, colnames(testres$tstat) !=
"Untagged_150_tube_vs_compl_Untagged_150_tube.t"]
tstat <- testres$tstat
colnames(tstat) <- .getProteinNameFromComparison(colnames(tstat))
TFs <- colnames(tstat)
## set NAs to minimum value
tstat[is.na(tstat)] <- min(na.omit(tstat))
## set all values below 0 to 0
tstat[tstat < 0] = 0
## convert it back to a matrix
heatmap_data_mean <- as.matrix(tstat)
## cluster the rows (to avoid slow clustering in LinkedCharts)
dist_mat1 <- dist(heatmap_data_mean, method = "euclidian")
hclust_avg1 <- hclust(dist_mat1, method = "average")
heatmap_data_mean_clustered <- heatmap_data_mean[hclust_avg1$order, ]
## calculate row max
heatmapRowMax <- apply(heatmap_data_mean_clustered, 1, max)
## add Pombase IDs and names and uniprot IDs
PomIds2 <- left_join(data.frame(row.names(heatmap_data_mean_clustered)),
idmap |> dplyr::mutate(unique_einprot_id = .capitalize(unique_einprot_id)),
by = c("row.names.heatmap_data_mean_clustered." =
"unique_einprot_id"))
PomIds1 <- ifelse(is.na(PomIds2$gene_stable_id),
PomIds2$row.names.heatmap_data_mean_clustered.,
PomIds2$gene_stable_id)
UniprotIDs <- PomIds2$xref
## convert the gene names to protein names
colnames(heatmap_data_mean_clustered) <-
.capitalize(colnames(heatmap_data_mean_clustered))
rownames(heatmap_data_mean_clustered) <-
.capitalize(rownames(heatmap_data_mean_clustered))
## save heatmap and related data as JSON
heatmap_data_mean2 <- toJSON(list(
heatmapMatrix = heatmap_data_mean_clustered,
TFs = colnames(heatmap_data_mean_clustered),
proteins = rownames(heatmap_data_mean_clustered),
MaxTSTAT = heatmapRowMax,
PomIds = PomIds1,
UniprotIDs = UniprotIDs),
pretty = TRUE)
write(heatmap_data_mean2, file = outfileHeatmap)
## Volcano plots
## -------------------------------------------------------------------------
## extract comparisons (include only plate format for Atf1 and untagged)
volcano_list <- metadata(sce)$testres$tests
volcano_list <- volcano_list[sub("\\.t", "", colnames(testres$tstat))]
TFIDs <- names(volcano_list)
TFnames <- .getProteinNameFromComparison(TFIDs)
# TFnames <- ifelse(grepl("^SP", TFnames), TFnames, stringr::str_to_title(TFnames))
## using sce object
vp_list <- lapply(
structure(colnames(heatmap_data_mean_clustered),
names = colnames(heatmap_data_mean_clustered)),
function(selTF) {
VP <- volcano_list[[which(TFnames == selTF)]][, c("pid", "logFC", "adj.P.Val",
"mlog10p", "showInVolcano")]
VP$DE <- ifelse(VP$showInVolcano, "significant", "not significant")
## replace NAs with 0
VP <- replace_na(data = VP,
replace = list(pid = "unknown",
logFC = 0,
adj.P.Val = 1,
mlog10p = 0,
DE = "not significant"))
## add to list
VP[, c("pid", "logFC", "adj.P.Val", "mlog10p", "DE")]
})
## convert the gene names to protein names
for (i in seq_along(vp_list)) {
rownames(vp_list[[i]]) <- .capitalize(rownames(vp_list[[i]]))
vp_list[[i]]$pid <- rownames(vp_list[[i]])
}
## save volcano plot data as JSON
vp_list2 <- toJSON(vp_list, dataframe = "columns", pretty = TRUE)
write(vp_list2, file = outfileVolcano)
}
.prepDataForTFexplorer(sce = sce150,
outfileHeatmap = "data/TSTAT_matrix4heatmap_150mM.json",
outfileVolcano = "data/volcanoData_150mM.json",
idmap = idmap)
.prepDataForTFexplorer(sce = sce500,
outfileHeatmap = "data/TSTAT_matrix4heatmap_500mM.json",
outfileVolcano = "data/volcanoData_500mM.json",
idmap = idmap)Facts and figures
Identification of baits in IP-MS
table(baitclass$class)##
## 150and500 150only nobait
## 40 43 7High stringency
Here we use an adjusted p-value cutoff of 0.0002 and a log2 fold change cutoff of 1.
makeInteractorSummary <- function(sce, baitdata, adjp_cutoff, logfc_cutoff, idmap,
baits_to_keep, tbl_title) {
ni_tc <- .getTestCols(sce = sce, adjp_cutoff = adjp_cutoff,
logfc_cutoff = logfc_cutoff)
ni_idx <- .getOrigBaitNameFromComparison(colnames(ni_tc$tstat), idmap = idmap) %in%
baits_to_keep
ni_int <- ni_tc$interactor[, ni_idx]
ni_baits <- cbind(match(.getProteinNameFromComparison(colnames(ni_int)),
rownames(ni_int)),
seq_len(ncol(ni_int)))
ni_int[ni_baits] <- NA
## Total number of interactions
table(ni_int > 0)
## Number of TFs with interactions
table(colSums(ni_int, na.rm = TRUE) > 0)
## Total number of TFs
ncol(ni_int)
## Total number of TF-TF interactions (no self-interactions)
sum(baitdata$nbrIntMat$nbrTFs_wobait)
## Number of TFs involved in these TF-TF interactions
tmp <- baitdata$nbrIntMat |>
filter(nbrTFs_wobait > 0)
length(unique(c(tmp$bait, unlist(strsplit(tmp$TFs, ",")))))
## Bilateral interactions
nbr_bilateral <- tmp |>
select(bait, TFs) |>
separate_longer_delim(TFs, delim = ",") |>
mutate(o1 = paste(bait, TFs, sep = "-"),
o2 = paste(TFs, bait, sep = "-")) |>
filter(o1 %in% o2) |>
nrow()
data.frame(`Total number of TFs for which bait was pulled down` = ncol(ni_int),
`Total number of interactions (excluding bait itself)` = length(which(ni_int > 0)),
`Number of TFs with interactions` = length(which(colSums(ni_int, na.rm = TRUE) > 0)),
`Number of TFs without observed interactions` = length(which(colSums(ni_int, na.rm = TRUE) == 0)),
`Total number of TF-TF interactions (excluding bait itself)` = sum(baitdata$nbrIntMat$nbrTFs_wobait),
`Number of TFs involved in these TF-TF interactions` = length(unique(c(tmp$bait, unlist(strsplit(tmp$TFs, ","))))),
`Number of bilateral TF-TF interaction` = nbr_bilateral,
check.names = FALSE) |>
pivot_longer(cols = everything(), names_to = "Aspect", values_to = "Value") |>
kbl(col.names = NULL, caption = tbl_title) |>
kable_styling()
}
## 150mM
## --------------------------------------------------------------------------
makeInteractorSummary(sce = sce150, baitdata = baitdata_150, adjp_cutoff = adjpThreshold,
logfc_cutoff = log2fcThreshold, idmap = idmap,
baits_to_keep = baitclass$Gene_name[baitclass$class != "nobait"],
tbl_title = "150 mM, high stringency")| Total number of TFs for which bait was pulled down | 83 |
| Total number of interactions (excluding bait itself) | 352 |
| Number of TFs with interactions | 43 |
| Number of TFs without observed interactions | 40 |
| Total number of TF-TF interactions (excluding bait itself) | 21 |
| Number of TFs involved in these TF-TF interactions | 17 |
| Number of bilateral TF-TF interaction | 16 |
## 500mM
## --------------------------------------------------------------------------
makeInteractorSummary(sce = sce500, baitdata = baitdata_500, adjp_cutoff = adjpThreshold,
logfc_cutoff = log2fcThreshold, idmap = idmap,
baits_to_keep = baitclass$Gene_name[baitclass$class != "nobait"],
tbl_title = "500 mM, high stringency")| Total number of TFs for which bait was pulled down | 40 |
| Total number of interactions (excluding bait itself) | 110 |
| Number of TFs with interactions | 24 |
| Number of TFs without observed interactions | 16 |
| Total number of TF-TF interactions (excluding bait itself) | 15 |
| Number of TFs involved in these TF-TF interactions | 14 |
| Number of bilateral TF-TF interaction | 12 |
Low stringency
Here we use an adjusted p-value cutoff of 0.01 and a log2 fold change cutoff of 1.
## 150mM
## --------------------------------------------------------------------------
makeInteractorSummary(sce = sce150, baitdata = baitdata_150_loose, adjp_cutoff = adjpThresholdLoose,
logfc_cutoff = log2fcThreshold, idmap = idmap,
baits_to_keep = baitclass$Gene_name[baitclass$class != "nobait"],
tbl_title = "150 mM, low stringency")| Total number of TFs for which bait was pulled down | 83 |
| Total number of interactions (excluding bait itself) | 881 |
| Number of TFs with interactions | 60 |
| Number of TFs without observed interactions | 23 |
| Total number of TF-TF interactions (excluding bait itself) | 39 |
| Number of TFs involved in these TF-TF interactions | 30 |
| Number of bilateral TF-TF interaction | 18 |
## 500mM
## --------------------------------------------------------------------------
makeInteractorSummary(sce = sce500, baitdata = baitdata_500_loose, adjp_cutoff = adjpThresholdLoose,
logfc_cutoff = log2fcThreshold, idmap = idmap,
baits_to_keep = baitclass$Gene_name[baitclass$class != "nobait"],
tbl_title = "500 mM, low stringency")| Total number of TFs for which bait was pulled down | 40 |
| Total number of interactions (excluding bait itself) | 202 |
| Number of TFs with interactions | 29 |
| Number of TFs without observed interactions | 11 |
| Total number of TF-TF interactions (excluding bait itself) | 17 |
| Number of TFs involved in these TF-TF interactions | 14 |
| Number of bilateral TF-TF interaction | 14 |
Figure 4
Settings
colLabel_150 <- expand.grid(c("Untagged", "Ace2", "Rst2", "Pap1"),
c("tube", "plate")) |>
tidyr::unite(col = "col", Var1, Var2) |>
dplyr::pull(col)
colLabel_500 <- expand.grid(c("Untagged", "Atf1", "Moc3", "Ntu1",
"Pcr1", "Ntu2"),
c("tube", "plate")) |>
tidyr::unite(col = "col", Var1, Var2) |>
dplyr::pull(col)
rowLabel <- NULL
complexLabel <- c("NuA4", "SAGA")Summary heatmaps using results from the statistical tests (150 mM NaCl)
Overall heatmap
hm1a <- Heatmap(hmdata_150$tstat,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = makeHeatmapCol(hmdata_150$tstat, stringency = "high"),
border = TRUE,
border_gp = gpar(lwd = 0.5),
column_split = hmdata_150$colsplit,
name = "Mod.\nt-stat.",
na_col = binary_heatmap_colors["FALSE"],
column_names_gp = gpar(fontsize = fs),
row_names_gp = gpar(fontsize = fs),
row_title = "Proteins significantly enriched in at least one experiment",
row_title_gp = gpar(fontsize = fl),
show_row_names = FALSE, show_column_names = FALSE,
use_raster = FALSE, show_heatmap_legend = TRUE,
right_annotation = makeComplexAnnotation(
hmdata_150$tstat, complexes[complexLabel],
idmap = idmap),
bottom_annotation = makeColLabels(hmdata_150$tstat, colLabel_150, fl),
width = unit(w_ipmsHm, "mm"), # heatmap body width
heatmap_height = unit(h_ipmsHm, "mm"), # whole heatmap height
column_title_gp = gpar(fontsize = ft),
heatmap_legend_param = list(title_gp = gpar(fontsize = fl),
border = "gray10",
border_gp = gpar(lwd = 0.5)))
set.seed(42L)
hm_150 <- grid.grabExpr(
hm1a <- draw(hm1a, merge_legend = TRUE,
heatmap_legend_side = "right",
annotation_legend_side = "right"),
width = w_ipmsHm, height = h_ipmsHm
)
plot_grid(hm_150)
Bait-vs-bait heatmap
nbrInteractors <- HeatmapAnnotation(
which = "column",
`Number of\ninteractors` =
anno_barplot(x = rowSums(as.matrix(baitdata_150$nbrIntMat[colnames(baitdata_150$mat),
c("nbrTFs_samefam",
"nbrTFs_difffam")])),
gp = gpar(fill = na_color,
col = na_color),
border = FALSE, bar_width = 1.0,
axis_param = list(gp = gpar(fontsize = fs)),
width = unit(10, "mm")),
annotation_name_side = "left",
annotation_name_rot = 0,
annotation_name_gp = gpar(fontsize = fl))
hm1b <- Heatmap(baitdata_150$mat,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = makeHeatmapCol(baitdata_150$mat, stringency = "high"),
border = TRUE,
border_gp = gpar(lwd = 0.5),
name = "Mod.\nt-stat.",
na_col = binary_heatmap_colors["FALSE"],
column_title = paste(c("150 mM NaCl", rep(" ", 45)), collapse = ""),
row_title = "Transcription factors",
column_title_gp = gpar(fontsize = fl + 3, fontface = "bold"),
row_title_gp = gpar(fontsize = fl),
column_names_gp = gpar(fontsize = fs),
row_names_gp = gpar(fontsize = fs),
show_row_names = FALSE, show_column_names = FALSE,
use_raster = FALSE, show_heatmap_legend = TRUE,
left_annotation = makeComplexAndDBDAnnotation(
baitdata_150$mat,
complexes[c("MBF transcription complex",
"CCAAT-binding factor complex",
"Atf1-Pcr1")],
dbd = dbd,
idmap = idmap),
top_annotation = nbrInteractors,
width = unit(w_ipmsHm, "mm") / 1.05, # heatmap body width
heatmap_height = unit(h_ipmsHm / 2.2, "mm"),
heatmap_legend_param = list(title_gp = gpar(fontsize = fl),
border = "gray10",
border_gp = gpar(lwd = 0.5)))
set.seed(42L)
hm_150_bait <- grid.grabExpr(
hm1b <- draw(hm1b, merge_legend = TRUE,
heatmap_legend_side = "right",
annotation_legend_side = "right"),
width = w_ipmsHm / 1.05, height = h_ipmsHm / 2.2
)
plot_grid(hm_150_bait)
## Separate legends
use_complexes <- c("MBF transcription complex",
"CCAAT-binding factor complex",
"Atf1-Pcr1")
lgd_complex <- Legend(labels = use_complexes,
legend_gp = gpar(fill = vapply(complexes[use_complexes],
function(x) x$color, NA_character_)),
title = "Complex", nrow = 2)
dbd_colors <- c(`KilA-N` = main_colors[5],
`Zn(II)2Cys6` = main_colors[4],
bZIP = main_colors[3],
`Other (combined)` = bg_color)
lgd_dbd <- Legend(labels = names(dbd_colors),
legend_gp = gpar(fill = dbd_colors),
title = "DBD family", nrow = 2)
legend_bait_1 <- grid.grabExpr(
hm1b <- ComplexHeatmap::draw(lgd_complex),
width = w_ipmsHm, height = h_ipmsHm / 8
)
legend_bait_2 <- grid.grabExpr(
hm1b <- ComplexHeatmap::draw(lgd_dbd),
width = w_ipmsHm, height = h_ipmsHm / 8
)Bait-vs-bait, less stringent adjusted p-value threshold
nbrInteractors <- HeatmapAnnotation(
which = "column",
`Number of\ninteractors` =
anno_barplot(x = rowSums(as.matrix(baitdata_150_loose$nbrIntMat[colnames(baitdata_150_loose$mat),
c("nbrTFs_samefam",
"nbrTFs_difffam")])),
gp = gpar(fill = na_color,
col = na_color),
border = FALSE, bar_width = 1.0,
axis_param = list(gp = gpar(fontsize = fs)),
width = unit(10, "mm")),
annotation_name_side = "left",
annotation_name_rot = 0,
annotation_name_gp = gpar(fontsize = fl))
hm1b_loose <- Heatmap(baitdata_150_loose$mat,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = makeHeatmapCol(baitdata_150_loose$mat, stringency = "low"),
border = TRUE,
border_gp = gpar(lwd = 0.5),
name = "Mod.\nt-stat.",
na_col = binary_heatmap_colors["FALSE"],
column_title = paste(c("150 mM NaCl", rep(" ", 45)), collapse = ""),
row_title = "Transcription factors",
column_title_gp = gpar(fontsize = fl + 3, fontface = "bold"),
row_title_gp = gpar(fontsize = fl),
column_names_gp = gpar(fontsize = fs),
row_names_gp = gpar(fontsize = fs),
show_row_names = FALSE, show_column_names = FALSE,
use_raster = FALSE, show_heatmap_legend = TRUE,
left_annotation = makeComplexAndDBDAnnotation(
baitdata_150_loose$mat,
complexes[c("MBF transcription complex",
"CCAAT-binding factor complex",
"Atf1-Pcr1")],
dbd = dbd,
idmap = idmap),
top_annotation = nbrInteractors,
width = unit(w_ipmsHm, "mm") / 1.05, # heatmap body width
heatmap_height = unit(h_ipmsHm / 2.2, "mm"),
heatmap_legend_param = list(title_gp = gpar(fontsize = fl),
border = "gray10",
border_gp = gpar(lwd = 0.5)))
hm_150_0.05_bait <- grid.grabExpr(
hm1b <- draw(hm1b_loose, merge_legend = TRUE,
heatmap_legend_side = "right",
annotation_legend_side = "right"),
width = w_ipmsHm / 1.05, height = h_ipmsHm / 2.2
)
plot_grid(hm_150_0.05_bait)
SAGA/NuA4 heatmap
tc150 <- .getTestCols(sce150, adjp_cutoff = 0.05,
logfc_cutoff = log2fcThreshold)
keep_proteins <- .capitalize(c(complexes$SAGA$gene_names, complexes$NuA4$gene_names))
keep_exp <- which(.getProteinNameFromComparison(colnames(tc150$adjp)) %in%
c("Ace2", "Pap1", "Rst2"))
binary_mat <- tc150$adjp[keep_proteins, keep_exp] < 0.05 &
tc150$logfc[keep_proteins, keep_exp] > 1
binary_mat[is.na(binary_mat)] <- FALSE
colnames(binary_mat) <- .getProteinNameFromComparison(colnames(binary_mat))
ordr <- order(rownames(binary_mat) %in% .capitalize(complexes$SAGA$gene_names),
binary_mat[, "Ace2"], binary_mat[, "Pap1"], binary_mat[, "Rst2"],
decreasing = TRUE)
binary_mat <- binary_mat[ordr, ]
binary_annot <- data.frame(NuA4 = rownames(binary_mat) %in%
.capitalize(complexes$NuA4$gene_names),
SAGA = rownames(binary_mat) %in%
.capitalize(complexes$SAGA$gene_names),
row.names = rownames(binary_mat))
mode(binary_mat) <- "character"
ht_opt(HEATMAP_LEGEND_PADDING = unit(0, "mm"))
hm1saga_nua4 <- Heatmap(
binary_mat,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = binary_heatmap_colors,
border = TRUE,
border_gp = gpar(lwd = 0.5),
name = " ",
na_col = binary_heatmap_colors["FALSE"],
column_names_gp = gpar(fontsize = ft),
row_names_gp = gpar(fontsize = fl),
show_row_names = TRUE, show_column_names = TRUE,
use_raster = FALSE, show_heatmap_legend = TRUE,
right_annotation = rowAnnotation(
df = binary_annot,
col = list(SAGA = c(`TRUE` = complex_colors[["SAGA"]],
`FALSE` = binary_heatmap_colors[["FALSE"]]),
NuA4 = c(`TRUE` = complex_colors[["NuA4"]],
`FALSE` = binary_heatmap_colors[["FALSE"]])),
show_legend = FALSE),
heatmap_height = unit(h_ipmsHm * 0.945, "mm"), # whole heatmap height
column_title_gp = gpar(fontsize = ft),
heatmap_legend_param = list(title_gp = gpar(fontsize = ft),
at = c("FALSE", "TRUE"),
labels = c("", "adj. p < 0.05 &\nlog2FC > 1"),
legend_gap = unit(0, "mm"))
)
hm_150_saga_nua4 <- grid.grabExpr(
hm1sn <- draw(hm1saga_nua4, merge_legend = TRUE,
heatmap_legend_side = "bottom",
annotation_legend_side = "right",
align_heatmap_legend = "global_center",
padding = unit(c(3, 1, 1, 1), "mm")),
width = w_ipmsHm, height = h_ipmsHm
)
plot_grid(hm_150_saga_nua4)
ht_opt(RESET = TRUE)RNA / DNA binding proteins
rnabinding <- read.delim(params$rnabinders)
keep <- rownames(hmdata_150$tstat) %in% c(rnabinding$PomBaseID, .capitalize(rnabinding$Gene_name))
sum(keep)## [1] 29hm1rna <- Heatmap(hmdata_150$tstat[keep, ],
cluster_rows = FALSE,
cluster_columns = TRUE,
col = makeHeatmapCol(hmdata_150$tstat, stringency = "high"),
border = TRUE,
border_gp = gpar(lwd = 0.5),
name = "Mod.\nt-stat.",
na_col = binary_heatmap_colors["FALSE"],
column_names_gp = gpar(fontsize = fs),
row_names_gp = gpar(fontsize = fs),
show_row_names = TRUE, show_column_names = FALSE,
use_raster = FALSE, show_heatmap_legend = TRUE,
# bottom_annotation = makeColLabels(tstat150[keep, ], colLabel_150, fl),
#left_annotation = makeRowLabels(tstat150[keep, ], rowLabel, fl),
width = unit(w_ipmsHm, "mm"), # heatmap body width
heatmap_height = unit(h_ipmsHm/2, "mm"), # whole heatmap height
heatmap_legend_param = list(title_gp = gpar(fontsize = fl),
border = "gray10",
border_gp = gpar(lwd = 0.5)))
set.seed(42L)
hm_150_rna <- grid.grabExpr(
hm1b <- draw(hm1rna, merge_legend = TRUE,
heatmap_legend_side = "right",
annotation_legend_side = "right"),
width = w_ipmsHm, height = h_ipmsHm/2
)
plot_grid(hm_150_rna)
dnabinding <- read.delim(params$dnabinders)
keep <- rownames(hmdata_150$tstat) %in% c(dnabinding$PomBaseID, .capitalize(dnabinding$Gene_name))
sum(keep)## [1] 127hm1dna <- Heatmap(hmdata_150$tstat[keep, ],
cluster_rows = FALSE,
cluster_columns = TRUE,
col = makeHeatmapCol(hmdata_150$tstat, stringency = "high"),
border = TRUE,
border_gp = gpar(lwd = 0.5),
name = "Mod.\nt-stat.",
na_col = binary_heatmap_colors["FALSE"],
column_names_gp = gpar(fontsize = fs),
row_names_gp = gpar(fontsize = fs),
show_row_names = TRUE, show_column_names = FALSE,
use_raster = FALSE, show_heatmap_legend = TRUE,
# bottom_annotation = makeColLabels(tstat150[keep, ], colLabel_150, fl),
#left_annotation = makeRowLabels(tstat150[keep, ], rowLabel, fl),
width = unit(w_ipmsHm, "mm"), # heatmap body width
heatmap_height = unit(1.5 * h_ipmsHm, "mm"), # whole heatmap height
heatmap_legend_param = list(title_gp = gpar(fontsize = fl),
border = "gray10",
border_gp = gpar(lwd = 0.5)))
set.seed(42L)
hm_150_dna <- grid.grabExpr(
hm1b2 <- draw(hm1dna, merge_legend = TRUE,
heatmap_legend_side = "right",
annotation_legend_side = "right"),
width = w_ipmsHm, height = 1.5 * h_ipmsHm
)
plot_grid(hm_150_dna)
Summary heatmaps using results from the statistical tests (500 mM NaCl)
Overall heatmap
hm2a <- Heatmap(hmdata_500$tstat,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = makeHeatmapCol(hmdata_500$tstat, stringency = "high"),
border = TRUE,
border_gp = gpar(lwd = 0.5),
column_split = hmdata_500$colsplit,
name = "Mod.\nt-stat.",
na_col = binary_heatmap_colors["FALSE"],
column_names_gp = gpar(fontsize = fs),
row_names_gp = gpar(fontsize = fs),
row_title = "Proteins significantly enriched in at least one experiment",
row_title_gp = gpar(fontsize = fl),
show_row_names = FALSE, show_column_names = FALSE,
use_raster = FALSE, show_heatmap_legend = TRUE,
right_annotation = makeComplexAnnotation(hmdata_500$tstat, complexes[complexLabel],
idmap = idmap, show_legend = FALSE),
bottom_annotation = makeColLabels(hmdata_500$tstat, colLabel_500, fl),
width = unit(w_ipmsHm, "mm"), # heatmap body width
heatmap_height = unit(h_ipmsHm, "mm"), # whole heatmap height
column_title_gp = gpar(fontsize = ft),
heatmap_legend_param = list(title_gp = gpar(fontsize = fl),
border = "gray10",
border_gp = gpar(lwd = 0.5)))
hm_500 <- grid.grabExpr(
hm2a <- draw(hm2a, merge_legend = TRUE,
heatmap_legend_side = "right",
annotation_legend_side = "right"),
width = w_ipmsHm, height = h_ipmsHm
)
plot_grid(hm_500)
Bait-vs-bait heatmap
nbrInteractors <- HeatmapAnnotation(
which = "column",
`Number of\ninteractors` =
anno_barplot(x = rowSums(as.matrix(baitdata_500$nbrIntMat[colnames(baitdata_500$mat),
c("nbrTFs_samefam",
"nbrTFs_difffam")])),
gp = gpar(fill = na_color,
col = na_color),
border = FALSE, bar_width = 1.0,
axis_param = list(gp = gpar(fontsize = fs)),
width = unit(10, "mm")),
annotation_name_side = "left",
annotation_name_rot = 0,
annotation_name_gp = gpar(fontsize = fl))
hm1b2 <- Heatmap(baitdata_500$mat,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = makeHeatmapCol(baitdata_500$mat, stringency = "high"),
border = TRUE,
border_gp = gpar(lwd = 0.5),
name = "Mod.\nt-stat.",
na_col = binary_heatmap_colors["FALSE"],
column_title = paste(c("500 mM NaCl", rep(" ", 45)), collapse = ""),
row_title = "Transcription factors",
column_title_gp = gpar(fontsize = fl + 3, fontface = "bold"),
row_title_gp = gpar(fontsize = fl),
column_names_gp = gpar(fontsize = fs),
row_names_gp = gpar(fontsize = fs),
show_row_names = FALSE, show_column_names = FALSE,
use_raster = FALSE, show_heatmap_legend = TRUE,
left_annotation = makeComplexAndDBDAnnotation(
baitdata_500$mat,
complexes[c("MBF transcription complex",
"CCAAT-binding factor complex", "Atf1-Pcr1")],
dbd = dbd, idmap = idmap),
top_annotation = nbrInteractors,
width = unit(w_ipmsHm, "mm") / 1.05, # heatmap body width
heatmap_height = unit(h_ipmsHm / 2.2, "mm"), # whole heatmap height
heatmap_legend_param = list(title_gp = gpar(fontsize = fl),
border = "gray10",
border_gp = gpar(lwd = 0.5)))
hm_500_bait <- grid.grabExpr(
hm1b2 <- draw(hm1b2, merge_legend = TRUE,
heatmap_legend_side = "right",
annotation_legend_side = "right"),
width = w_ipmsHm / 1.05, height = h_ipmsHm / 2.2
)
plot_grid(hm_500_bait)
Compare t-statistics from low-salt and high-salt conditions
The figures below show scatter plots of moderated t-statistics for selected bait pull-downs compared to their complements, in the high-salt vs low-salt conditions.
## Extract results from statistical tests
rdlow <- rowData(sce150)
rdhigh <- rowData(sce500)
## Get t-statistics
contrasts_to_compare <- grep("\\.t$", colnames(rdhigh), value = TRUE)
contrasts_low <- gsub("_compl", "_highsaltcompl",
gsub("_500_", "_150_", contrasts_to_compare))
idx <- which(!contrasts_low %in% colnames(rdlow))
contrasts_low[idx] <- gsub("_plate", "_tube", contrasts_low[idx])
stopifnot(all(contrasts_low %in% colnames(rdlow)))
## t-statistics, logFCs and adjusted p-values, low salt
low <- as.data.frame(rdlow[, c("einprotId", contrasts_low)]) |>
tidyr::pivot_longer(names_to = "contrast_low", values_to = "tstat_low", -einprotId) |>
dplyr::mutate(contrast_low = sub("\\.t$", "", contrast_low)) |>
dplyr::left_join(
as.data.frame(rdlow[, c("einprotId", sub("\\.t$", ".adj.P.Val",
contrasts_low))]) |>
tidyr::pivot_longer(names_to = "contrast_low",
values_to = "adjp_low", -einprotId) |>
dplyr::mutate(contrast_low = sub("\\.adj.P.Val$", "", contrast_low))) |>
dplyr::left_join(
as.data.frame(rdlow[, c("einprotId", sub("\\.t$", ".logFC",
contrasts_low))]) |>
tidyr::pivot_longer(names_to = "contrast_low",
values_to = "logfc_low", -einprotId) |>
dplyr::mutate(contrast_low = sub("\\.logFC$", "", contrast_low)))
## t-statistics, logFCs and adjusted p-values, high salt
high <- as.data.frame(rdhigh[, c("einprotId", contrasts_to_compare)]) |>
tidyr::pivot_longer(names_to = "contrast_high",
values_to = "tstat_high", -einprotId) |>
dplyr::mutate(contrast_high = sub("\\.t$", "", contrast_high)) |>
dplyr::left_join(
as.data.frame(rdhigh[, c("einprotId", sub("\\.t$", ".adj.P.Val",
contrasts_to_compare))]) |>
tidyr::pivot_longer(names_to = "contrast_high",
values_to = "adjp_high", -einprotId) |>
dplyr::mutate(contrast_high = sub("\\.adj.P.Val$", "",
contrast_high))) |>
dplyr::left_join(
as.data.frame(rdhigh[, c("einprotId", sub("\\.t$", ".logFC",
contrasts_to_compare))]) |>
tidyr::pivot_longer(names_to = "contrast_high",
values_to = "logfc_high", -einprotId) |>
dplyr::mutate(contrast_high = sub("\\.logFC$", "",
contrast_high))) |>
dplyr::mutate(contrast_low = sub("\\.t$", "",
contrasts_low[match(paste0(contrast_high, ".t"),
contrasts_to_compare)]))
## Merge data.frames
lowhigh <- low |> dplyr::full_join(high) |>
dplyr::mutate(tstat_low = replace(tstat_low, is.na(tstat_low), 0),
tstat_high = replace(tstat_high, is.na(tstat_high), 0),
logfc_low = replace(logfc_low, is.na(logfc_low), 0),
logfc_high = replace(logfc_high, is.na(logfc_high), 0),
adjp_low = replace(adjp_low, is.na(adjp_low), 1),
adjp_high = replace(adjp_high, is.na(adjp_high), 1)) |>
dplyr::mutate(bait = .getProteinNameFromComparison(contrast_low)) |>
dplyr::mutate(splitvar = .getProteinNameFromComparison(contrast_low))
## Add dal2
rddal2 <- rowData(scedal2)
contrasts_low_dal2 <- c("Dal2_150_plate_vs_Untagged_150_plate.t")
low_dal2 <- as.data.frame(rddal2[, c("einprotId", contrasts_low_dal2)]) |>
tidyr::pivot_longer(names_to = "contrast_low",
values_to = "tstat_low", -einprotId) |>
dplyr::mutate(contrast_low = sub("\\.t$", "", contrast_low)) |>
dplyr::left_join(
as.data.frame(rddal2[, c("einprotId", sub("\\.t$", ".adj.P.Val",
contrasts_low_dal2))]) |>
tidyr::pivot_longer(names_to = "contrast_low",
values_to = "adjp_low", -einprotId) |>
dplyr::mutate(contrast_low = sub("\\.adj.P.Val$", "", contrast_low))) |>
dplyr::left_join(
as.data.frame(rddal2[, c("einprotId", sub("\\.t$", ".logFC",
contrasts_low_dal2))]) |>
tidyr::pivot_longer(names_to = "contrast_low",
values_to = "logfc_low", -einprotId) |>
dplyr::mutate(contrast_low = sub("\\.logFC$", "", contrast_low)))
contrasts_high_dal2 <- c("Dal2_500_plate_vs_Untagged_500_plate.t")
high_dal2 <- as.data.frame(rddal2[, c("einprotId", contrasts_high_dal2)]) |>
tidyr::pivot_longer(names_to = "contrast_high",
values_to = "tstat_high", -einprotId) |>
dplyr::mutate(contrast_high = sub("\\.t$", "", contrast_high)) |>
dplyr::left_join(
as.data.frame(rddal2[, c("einprotId", sub("\\.t$", ".adj.P.Val",
contrasts_high_dal2))]) |>
tidyr::pivot_longer(names_to = "contrast_high",
values_to = "adjp_high", -einprotId) |>
dplyr::mutate(contrast_high = sub("\\.adj.P.Val$", "", contrast_high))) |>
dplyr::left_join(
as.data.frame(rddal2[, c("einprotId", sub("\\.t$", ".logFC",
contrasts_high_dal2))]) |>
tidyr::pivot_longer(names_to = "contrast_high",
values_to = "logfc_high", -einprotId) |>
dplyr::mutate(contrast_high = sub("\\.logFC$", "", contrast_high)))
lowhigh_dal2 <- low_dal2 |> dplyr::full_join(high_dal2) |>
dplyr::mutate(tstat_low = replace(tstat_low, is.na(tstat_low), 0),
tstat_high = replace(tstat_high, is.na(tstat_high), 0),
logfc_low = replace(logfc_low, is.na(logfc_low), 0),
logfc_high = replace(logfc_high, is.na(logfc_high), 0),
adjp_low = replace(adjp_low, is.na(adjp_low), 1),
adjp_high = replace(adjp_high, is.na(adjp_high), 1)) |>
dplyr::mutate(bait = .getProteinNameFromComparison(sub("SPB4775_", "",
contrast_low))) |>
dplyr::mutate(splitvar = .getProteinNameFromComparison(sub("SPB4775_", "",
contrast_low)))
lowhigh <- lowhigh |>
dplyr::bind_rows(lowhigh_dal2)
## Subset 1
adjp150vs500 <- 1e-3
logfc150vs500 <- 1
lowhighsub1 <- lowhigh |>
dplyr::filter(splitvar %in% c("Ace2", "Rst2", "Pap1")) |>
dplyr::mutate(splitvar = factor(splitvar, levels = c("Ace2", "Rst2", "Pap1"))) |>
dplyr::mutate(significant = (adjp_high < adjp150vs500 & logfc_high > logfc150vs500)) |>
dplyr::mutate(complex = ifelse(einprotId %in% .capitalize(complexes$NuA4$gene_names),
"NuA4", ifelse(einprotId %in% .capitalize(complexes$SAGA$gene_names),
"SAGA", NA_character_)))
## Subset 2
lowhighsub2 <- lowhigh |>
dplyr::filter(splitvar %in% c("Atf1", "Pcr1", "Moc3",
"Dal2", "Ntu1", "Ntu2")) |>
dplyr::mutate(splitvar = factor(splitvar, levels = c("Atf1", "Pcr1", "Moc3",
"Dal2", "Ntu1", "Ntu2"))) |>
dplyr::mutate(significant = (adjp_high < adjp150vs500 & logfc_high > logfc150vs500))
(gg_retain1 <- ggplot(lowhighsub1,
aes(x = tstat_high, y = tstat_low)) +
geom_abline(data = lowhighsub1 |>
dplyr::filter(einprotId == bait),
aes(slope = tstat_low / tstat_high, intercept = 0),
linetype = "dashed", color = "grey50") +
geom_text_repel(data = lowhighsub1 |>
dplyr::mutate(
einprotId = ifelse(significant, einprotId, "")),
aes(label = einprotId), min.segment.length = 0,
max.overlaps = Inf) +
geom_point(aes(alpha = !is.na(complex), color = complex),
show.legend = c(alpha = FALSE, color = TRUE)) +
scale_color_manual(values = complex_colors, na.value = "black",
breaks = c("SAGA", "NuA4")) +
scale_alpha_manual(values = c(`TRUE` = 1, `FALSE` = 0.25)) +
geom_point(data = lowhighsub1 |> dplyr::filter(einprotId == bait) |>
dplyr::mutate(type = "Bait"),
mapping = aes(fill = type), color = main_colors[3],
shape = 21, alpha = 1, size = 2) +
scale_fill_manual(values = c(Bait = main_colors[3])) +
facet_wrap(~ splitvar, scales = "free", ncol = 1,
labeller = labeller(splitvar = structure(
paste0(unique(lowhighsub1$splitvar), " IP-MS"),
names = as.character(unique(lowhighsub1$splitvar))))) +
labs(x = "500 mM NaCl\nmod. t-stat.", y = "150 mM NaCl mod. t-stat.") +
theme_cowplot(ft) +
theme(legend.position = "bottom",
strip.background = element_rect(fill = NA, colour = NA),
strip.text = element_text(hjust = 0, face = "bold")) +
guides(color = guide_legend(nrow = 1, byrow = TRUE, title = ""),
fill = guide_legend(nrow = 1, byrow = TRUE, title = "")))
(gg_retain2 <- ggplot(lowhighsub2,
aes(x = tstat_high, y = tstat_low)) +
geom_abline(data = lowhighsub2 |>
dplyr::filter(einprotId == bait),
aes(slope = tstat_low / tstat_high, intercept = 0),
linetype = "dashed", color = "grey50") +
geom_text_repel(data = lowhighsub2 |>
dplyr::mutate(einprotId = ifelse(significant, einprotId, "")),
aes(label = einprotId), min.segment.length = 0,
max.overlaps = Inf) +
geom_point(aes(alpha = significant, color = significant, size = significant),
show.legend = c(alpha = FALSE, color = TRUE, size = FALSE)) +
scale_color_manual(
values = c(`TRUE` = lighten(main_colors[4], 0.2), `FALSE` = "black"),
labels = c(`TRUE` = paste0("adj. p < ", adjp150vs500,
" and logFC > ", logfc150vs500, " (500 mM NaCl)"),
`FALSE` = paste0("adj. p >= ", adjp150vs500,
" or logFC <= ", logfc150vs500, " (500 mM NaCl)"))) +
scale_size_manual(values = c(`TRUE` = 2, `FALSE` = 1)) +
scale_alpha_manual(values = c(`TRUE` = 1, `FALSE` = 0.25)) +
geom_point(data = lowhighsub2 |>
dplyr::filter(einprotId == bait) |>
dplyr::mutate(type = "Bait"),
mapping = aes(fill = type), color = main_colors[3],
shape = 21, alpha = 1, size = 2) +
scale_fill_manual(values = c(Bait = main_colors[3])) +
facet_wrap(~ splitvar, scales = "free", ncol = 2,
labeller = labeller(splitvar = structure(
paste0(unique(lowhighsub2$splitvar), " IP-MS"),
names = as.character(unique(lowhighsub2$splitvar))))) +
labs(x = "500 mM NaCl mod. t-stat.", y = "150 mM NaCl mod. t-stat.") +
theme_cowplot(ft) +
theme(legend.position = "bottom",
strip.background = element_rect(fill = NA, colour = NA),
strip.text = element_text(hjust = 0, face = "bold")) +
guides(color = guide_legend(nrow = 2, byrow = TRUE, title = ""),
fill = guide_legend(nrow = 1, byrow = TRUE, title = "")))
Supplementary plots - PCA
## Low salt
sce150_pca <- einprot::doPCA(
sce150[rowSums(!assay(sce150,
metadata(sce150)$aNames$assayImputIndic)) >= 2, ],
metadata(sce150)$aNames$assayImputed)$sce
df <- scuttle::makePerCellDF(sce150_pca, features = NULL, use.coldata = TRUE,
use.dimred = TRUE) |>
dplyr::mutate(Protocol = .capitalize(sub("_digest|_freeze", "", sub("SOP_", "", Protocol))))
pcvar150 <- attr(reducedDim(sce150_pca, "PCA_log2_LFQ.intensity"), "percentVar")
(ggpca1 <- ggplot(
df,
aes(x = .data[[paste0("PCA_", metadata(sce150)$aNames$assayImputed, ".1")]],
y = .data[[paste0("PCA_", metadata(sce150)$aNames$assayImputed, ".2")]],
color = Protocol)) +
geom_point(size = 3, alpha = 0.5) +
scale_color_manual(values = c(Plate = main_colors[1],
Tube = main_colors[5]), name = "") +
coord_fixed() +
labs(x = paste0("PC1 (", round(pcvar150[1], 1), "%)"),
y = paste0("PC2 (", round(pcvar150[2], 1), "%)"),
title = "150 mM NaCl") +
theme_cowplot() +
theme(legend.position = "bottom") +
guides(color = guide_legend(nrow = 1, byrow = TRUE)))
(ggpca2 <- ggplot(
df |>
dplyr::mutate(bait = ifelse(Gene_name == "atf1", "Atf1", "Other baits")) |>
dplyr::mutate(prot = .capitalize(prot)),
aes(x = .data[[paste0("PCA_", metadata(sce150)$aNames$assayImputed, ".1")]],
y = .data[[paste0("PCA_", metadata(sce150)$aNames$assayImputed, ".2")]],
color = bait, shape = prot)) +
geom_point(size = 3, alpha = 0.5) +
scale_color_manual(values = c(Atf1 = main_colors[3],
`Other baits` = na_color), name = "") +
scale_shape_manual(values = c(Plate = 16, Tube = 17), name = "") +
coord_fixed() +
labs(x = paste0("PC1 (", round(pcvar150[1], 1), "%)"),
y = paste0("PC2 (", round(pcvar150[2], 1), "%)"), title = "150 mM NaCl") +
theme_cowplot() +
theme(legend.position = "bottom", legend.box = "horizontal",
legend.margin = margin()))
sce500_pca <- einprot::doPCA(
sce500[rowSums(!assay(sce500, metadata(sce500)$aNames$assayImputIndic)) >= 2, ],
metadata(sce500)$aNames$assayImputed)$sce
df <- scuttle::makePerCellDF(sce500_pca, features = NULL, use.coldata = TRUE,
use.dimred = TRUE)
pcvar500 <- attr(reducedDim(sce500_pca, "PCA_log2_LFQ.intensity"), "percentVar")
(ggpca3 <- ggplot(
df |>
dplyr::mutate(bait = ifelse(Gene_name == "atf1", "Atf1", "Other baits")),
aes(x = .data[[paste0("PCA_", metadata(sce500)$aNames$assayImputed, ".1")]],
y = .data[[paste0("PCA_", metadata(sce500)$aNames$assayImputed, ".2")]],
color = bait)) +
geom_point(size = 3, alpha = 0.5) +
scale_color_manual(values = c(Atf1 = main_colors[3],
`Other baits` = na_color), name = "") +
coord_fixed() +
labs(x = paste0("PC1 (", round(pcvar500[1], 1), "%)"),
y = paste0("PC2 (", round(pcvar500[2], 1), "%)"),
title = "500 mM NaCl") +
theme_cowplot() +
theme(legend.position = "bottom"))
Supplementary plot - tube vs plate t-statistics for Atf1 pull-down
The low-salt experiments were done using two different protocols (tube or plate). The Atf1 pull-down was performed with both protocols, and the figure below compares the moderated t-statistics (each obtained by comparing to the complement consisting of all non-Atf1 pull-down samples obtained with the same protocol).
df <- as.data.frame(rowData(sce150)) |>
dplyr::select(einprotId, Atf1_150_tube_vs_compl_Atf1_150_tube.t,
Atf1_150_plate_vs_compl_Atf1_150_plate.t) |>
dplyr::mutate(col = ifelse(einprotId == "Atf1", "bait",
ifelse(einprotId == "Pcr1", "Pcr1", "other")))
(ggatf1 <- ggplot(df, aes(x = Atf1_150_tube_vs_compl_Atf1_150_tube.t,
y = Atf1_150_plate_vs_compl_Atf1_150_plate.t)) +
geom_abline(slope = 1, intercept = 0) +
geom_point(size = 3, aes(color = col, alpha = col)) +
scale_color_manual(values = c(bait = main_colors[3], pcr1 = main_colors[4],
other = "black")) +
scale_alpha_manual(values = c(bait = 1, pcr1 = 1, other = 0.5)) +
geom_text_repel(data = df |> filter(einprotId %in% c("Atf1", "Pcr1")),
aes(label = einprotId)) +
labs(x = "Mod. t-stat. Atf1 vs complement\n(tube, 150 mM NaCl)",
y = "Mod. t-stat. Atf1 vs complement\n(plate, 150 mM NaCl)") +
coord_fixed() +
theme_cowplot() +
theme(legend.position = "none"))## Warning: Removed 557 rows containing missing values or values outside the scale range (`geom_point()`).
Put together
cowplot::plot_grid(
hm_150,
plot_grid(gg_retain1, NULL, hm_150_saga_nua4, nrow = 1, rel_widths = c(1, 0.05, 0.75)),
hm_500,
gg_retain2,
nrow = 2,
labels = c("A", "C", "B", "D"),
align = "v",
axis = "b"
)
Supplementary Figure
cowplot::plot_grid(
wrap_plots(ggpca1, ggpca2, ggatf1, ggpca3, hm_150_bait, hm_500_bait, ncol = 2) +
plot_annotation(tag_levels = list(c("A", "", "", "", "B"))) &
theme(plot.tag = element_text(size = 14, face = "bold")),
cowplot::plot_grid(
NULL,
legend_bait_1,
legend_bait_2,
NULL,
nrow = 1,
rel_widths = c(1, 2, 2, 1)
),
ncol = 1,
rel_heights = c(12, 1)
)
Session info
Session info
sessionInfo()## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: FlexiBLAS OPENBLAS; LAPACK version 3.10.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Zurich
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.2.0 kableExtra_1.4.0 tidyr_1.3.1 stringr_1.5.1
## [5] DT_0.33 ggrepel_0.9.5 colorspace_2.1-0 ComplexHeatmap_2.18.0
## [9] jsonlite_1.8.8 forcats_1.0.0 dplyr_1.1.4 einprot_0.9.4
## [13] cowplot_1.1.3 ggplot2_3.5.0 SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
## [17] Biobase_2.62.0 GenomicRanges_1.54.1 GenomeInfoDb_1.38.8 IRanges_2.36.0
## [21] S4Vectors_0.40.2 BiocGenerics_0.48.1 MatrixGenerics_1.14.0 matrixStats_1.3.0
##
## loaded via a namespace (and not attached):
## [1] R.methodsS3_1.8.2 ggalt_0.4.0 poweRlaw_0.80.0 Biostrings_2.70.3
## [5] vctrs_0.6.5 digest_0.6.35 png_0.1-8 shape_1.4.6.1
## [9] pcaPP_2.0-4 iSEE_2.14.0 magick_2.8.3 MASS_7.3-60
## [13] reshape2_1.4.4 httpuv_1.6.15 foreach_1.5.2 withr_3.0.0
## [17] xfun_0.43 survival_3.5-8 memoise_2.0.1 hexbin_1.28.3
## [21] ggbeeswarm_0.7.2 gmm_1.8 systemfonts_1.0.6 zoo_1.8-12
## [25] GlobalOptions_0.1.2 gtools_3.9.5 R.oo_1.26.0 DEoptimR_1.1-3
## [29] GGally_2.2.1 KEGGREST_1.42.0 promises_1.3.0 httr_1.4.7
## [33] restfulr_0.0.15 hash_2.2.6.3 rstudioapi_0.16.0 shinyAce_0.4.2
## [37] rrcovNA_0.5-1 miniUI_0.1.1.1 generics_0.1.3 babelgene_22.9
## [41] zlibbioc_1.48.2 ScaledMatrix_1.10.0 GenomeInfoDbData_1.2.11 SparseArray_1.2.4
## [45] xtable_1.8-4 ade4_1.7-22 pracma_2.4.4 doParallel_1.0.17
## [49] evaluate_0.23 ExploreModelMatrix_1.14.0 S4Arrays_1.2.1 proDA_1.16.0
## [53] hms_1.1.3 irlba_2.3.5.1 iSEEhex_1.4.0 shinyWidgets_0.8.5
## [57] magrittr_2.0.3 readr_2.1.5 later_1.3.2 viridis_0.6.5
## [61] lattice_0.22-6 MsCoreUtils_1.14.1 genefilter_1.84.0 robustbase_0.99-2
## [65] XML_3.99-0.16.1 scuttle_1.12.0 pillar_1.9.0 nlme_3.1-164
## [69] iterators_1.0.14 caTools_1.18.2 compiler_4.3.2 beachmat_2.18.1
## [73] stringi_1.8.3 GenomicAlignments_1.38.2 tmvtnorm_1.6 plyr_1.8.9
## [77] msigdbr_7.5.1 crayon_1.5.2 abind_1.4-5 BiocIO_1.12.0
## [81] scater_1.30.1 chron_2.3-61 bit_4.0.5 sandwich_3.1-0
## [85] pcaMethods_1.94.0 codetools_0.2-19 BiocSingular_1.18.0 crosstalk_1.2.1
## [89] bslib_0.7.0 GetoptLong_1.0.5 plotly_4.10.4 mime_0.12
## [93] MultiAssayExperiment_1.28.0 splines_4.3.2 circlize_0.4.16 Rcpp_1.0.12
## [97] sparseMatrixStats_1.14.0 Rttf2pt1_1.3.12 knitr_1.46 blob_1.2.4
## [101] utf8_1.2.4 clue_0.3-65 seqLogo_1.68.0 AnnotationFilter_1.26.0
## [105] QFeatures_1.12.0 DelayedMatrixStats_1.24.0 tibble_3.2.1 sqldf_0.4-11
## [109] iSEEu_1.14.0 Matrix_1.6-5 statmod_1.5.0 tzdb_0.4.0
## [113] svglite_2.1.3 pkgconfig_2.0.3 tools_4.3.2 cachem_1.0.8
## [117] RSQLite_2.3.6 viridisLite_0.4.2 DBI_1.2.2 impute_1.76.0
## [121] fastmap_1.1.1 rmarkdown_2.26 scales_1.3.0 shinydashboard_0.7.2
## [125] Rsamtools_2.18.0 sass_0.4.9 ggstats_0.6.0 farver_2.1.1
## [129] mgcv_1.9-1 gsubfn_0.7 yaml_2.3.8 rtracklayer_1.62.0
## [133] cli_3.6.2 purrr_1.0.2 writexl_1.5.0 lifecycle_1.0.4
## [137] mvtnorm_1.2-4 rintrojs_0.3.4 BiocParallel_1.36.0 annotate_1.80.0
## [141] gtable_0.3.5 rjson_0.2.21 parallel_4.3.2 limma_3.58.1
## [145] colourpicker_1.3.0 TFBSTools_1.40.0 bitops_1.0-7 bit64_4.0.5
## [149] BiocNeighbors_1.20.2 proto_1.0.0 CNEr_1.38.0 jquerylib_0.1.4
## [153] highr_0.10 shinyjs_2.1.0 imputeLCMD_2.1 R.utils_2.12.3
## [157] rrcov_1.7-5 lazyeval_0.2.2 shiny_1.8.1.1 htmltools_0.5.8.1
## [161] proj4_1.0-14 GO.db_3.18.0 glue_1.7.0 STRINGdb_2.14.3
## [165] TFMPvalue_0.0.9 XVector_0.42.0 RCurl_1.98-1.14 ComplexUpset_1.3.3
## [169] mclust_6.1 BSgenome_1.70.2 motifStack_1.46.0 gridExtra_2.3
## [173] igraph_2.0.3 extrafontdb_1.0 R6_2.5.1 ggiraph_0.8.9
## [177] gplots_3.1.3.1 labeling_0.4.3 cluster_2.1.4 stringdist_0.9.12
## [181] DirichletMultinomial_1.44.0 DelayedArray_0.28.0 tidyselect_1.2.1 vipor_0.4.7
## [185] plotrix_3.8-4 ProtGenerics_1.34.0 maps_3.4.2 xml2_1.3.6
## [189] ash_1.0-15 AnnotationDbi_1.64.1 rsvd_1.0.5 munsell_0.5.1
## [193] KernSmooth_2.23-22 data.table_1.15.4 htmlwidgets_1.6.4 RColorBrewer_1.1-3
## [197] rlang_1.1.3 extrafont_0.19 uuid_1.2-0 norm_1.0-11.1
## [201] fansi_1.0.6 Cairo_1.6-2 beeswarm_0.4.0